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Abstract 



Kinetic equations governing time evolution of positions and momenta of atoms in extended systems are de- 
O ' rived using quantum-classical ensembles within the Non-Equilibrium Statistical Operator Method (NESOM). 

Ions are treated classically, while their electrons quantum mechanically; however, the statistical operator is not 
factorised in any way and no simplifying assumptions are made concerning the electronic subsystem. Using 
this method, we derive kinetic equations of motion for the classical degrees of freedom (atoms) which account 
fully for the interaction and energy exchange with the quantum variables (electrons). Our equations, alongside 
' the usual "Newtonian"-like terms normally associated with the Ehrenfest dynamics, contain additional terms, 

proportional to the atoms velocities, which can be associated with the "electronic friction". Possible ways of 
calculating the friction forces which are shown to be given via complicated non-equilibrium correlation func- 
tions, are discussed. In particular, we demonstrate that the correlation functions are directly related to the 
thermodynamic Matsubara Green's functions, and this relationship allows for the diagrammatic methods to be 
used in treating electron-electron interaction perturbatively when calculating the correlation functions. This 
work also generalises previous attempts, mostly based on model systems, of introducing the electronic friction 
into Molecular Dynamics equations of atoms. 

1 Introduction 

> . 

Classical Molecular Dynamics (MD) simulations pQ,[2] play an important role in modern condensed matter physics [5] 
giving direct access to a wide range of statistical properties of the systems under study. In ab initio MD simulations 
atoms which are treated classically follow in time Newtonian's equations of motion. The latter are solved numerically 
using atomic forces. In ab initio MD simulations the forces on atoms are calculated from the first principles by 
considering electrons (at each atomic configuration) entirely quantum mechanically, usually within the density 
functional theory (DFT) [4]. This approach is also sometimes called the mean-field approximation (MFA) [5]. 

Probably, the simplest quantum-mechanical justification of the MFA [6] is based on a factorisation of the density 
operator for the whole system into a product of individual operators for the nuclei and electrons, and then using this 
Ansatz in the quantum Liouville equation with subsequent replacement of the quantum bracket with the classical 
Poisson bracket for the classical degrees of freedom. Then, a classical trajectory is introduced by adopting a special 
Delta-function representation for the density operator of the classical subsystem. The important message here is 
that the ionic coordinates and momenta in the usual MD equations appear as statistical averages calculated at 
every time step. Thus, usual MD equations constitute the dynamical equations of motion (EoM) for averages 
as proposed originally a long time ago by Ehrenfest [Z]- Although more sophisticated approaches have also been 
developed (see e.g. [HI El ED] and references therein) in which quantum nature of slow variables (classical degrees 
of freedom) is taken into account to some extent, these methods are still very complicated. At the same time, 
classical consideration of nuclei can still be well justified for many problems [TT] . 

In the present paper we propose a general statistical mechanical consideration of a system consisting of slow 
and fast degrees of freedom assuming nuclei and electrons as a particular example. We derive EoM for slow degrees 
of freedom (nuclei) which interact and exchange energy with the fast degrees of freedom (electrons) . Contrary to 
conventional approaches (e.g. [9]) based on the Liouville equation which possesses the time-reversal symmetry and 
is thus intrinsically equilibrium [32] , our method is based on entirely non-equilibrium consideration within the Non- 
equilibrium Statistical Operator Method (NESOM) [12]. The treatment of classical and quantum degrees of freedom 
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is done within the method of mixed quantum-classical ensembles (MQCE) [13| fT4] in which the Liouville operator is 
represented as a sum of the classical Poisson bracket and a quantum- mechanical commutator acting on the statistical 
operator of the whole system. The latter depends on coordinates and momenta of classical degrees of freedom and, 
at the same time, is a super-operator acting on quantum operators due to quantum degrees of freedom. A formal 
derivation of this method of treating mixed quantum-classical systems based on a group-theoretical analysis was 
given in Refs. [HI [US [17] . Note that in this approach the statistical operator is not assumed to be in a factorised 
form with respect to slow and fast degrees of freedom. 

Physically, one would expect that if the "fast" electrons are in instantaneous equilibrium with the "slow" nuclei, 
then the former would follow the dynamics of the latter and the EoM would correspond exactly to the Ehrenfest 
dynamics when the nuclei move along a single trajectory (which depends on the initial conditions) while the 
electrons are in the ground state. However, in reality the electrons are quantum particles which impose fluctuating 
forces on the nuclei. Also, at a given temperature T, the electrons are not isolated from the nuclei and the heat 
bath surrounding the system; they should get enough energy to occupy an ensemble of ground and excited states 
corrresponding to this T and will require some time to equilibrate after the nuclei displaced from their current 
positions. This make us think that the motion of the nuclei cannot be considered as following a single trajectory; 
instead, one can only consider the motion of the nuclei statistically, "on average". Moreover, the EoM for the nuclei 
average momenta would also deviate from the Ehrenfest dynamics: firstly, the average forces acting on the nuclei 
are expected to contain friction-like terms reflecting the possibility of the energy exchange, and, secondly, there 
should be a "conservative" force acting on the nuclei due to electrons occupying an ensemble of states. In this paper 
we develop a general formalism that leads to this kind of description. 

We show that the EoM for nuclei corresponds to the Ehrenfest dynamics with additional terms. The latter are 
related to rather complicated non-equilibrium correlation functions, and we provide a way of deriving these terms 
systematically. In the first order approximation, our additional terms are shown to be exactly proportional to atoms 
momenta and can thus be interpreted as friction forces. These forces have long been known in the literature as 
"electronic friction" (see review [6j and references therein), but they were either introduced semi-empirically [18j, 
as Langevin forces [H] or due to energy losses in particular model systems [fH HU [2Ql El] ■ In this paper we give 
a general derivation and justification of these kind of terms. The method presented here is a generalisation of 
our previous treatment |13l [T4] of a classical tip of Atomic Force Microscopy interacting with quantum surface 
vibrations. 

The plan of the paper is as follows. In the next Section we shall introduce main concepts of the NESOM and 
MQCE to set up the necessary definitions and notations. In Section 3 our main formalism is given and the EoM for 
the nuclei are derived to the first order, and we explain which quantities are used for building up this approximation. 
We also discuss how this procedure can be extended systematically to include terms up to any order. Although we 
do not consider any specific model in which the non-equilibrium correlation functions could be calculated due to 
an enormous complexity of those, a general discussion on how this could be done, at least in principle, will also be 
given. In particular, their connection to the Matsubara Green's functions ESIEHJ is discussed in Appendix. Finally, 
in Section 4, main conclusions are drawn. 



2 NESOM and MQCE 

In NESOM statistical mechanics of a system is described in general by a statistical operator p(t) which satisfies the 
Liouville equation with broken time-reversal symmetry |12| : 

dp 

— + iLp = -S (p - p re l) (1) 

where L is the Liouville operator and p re i (t) is the so-called relevant distribution corresponding to local equilibrium 
in the system. The right hand side of Eq. (TJ serves to break the time-reversal symmetry inherent to the usual 
Liouville equation, ^| + iLp = 0, in which this term is missing, and this guarantees that the retarded solution 
of the Liouville equation is chosen corresponding to the physically acceptable non-equilibrium behaviour of the 
system. The limit e — > (following the thermodynamic limit) is taken after the calculation of necessary averages 
with respect to p(t). 

The key quantity in this formalism is the relevant distribution, p re i{t), which is constructed from a set of 
relevant statistical variables X n , and follows from the principle of maximum of the information entropy. The 
entropy is maximum subject to the so-called self-consistency conditions stating that the statistical averages of 
the variables X n calculated with the relevant distribution, (X n ) rel = Tr (p re i(t)X n ), at time t are always equal 
to the true statistical averages, (X n ) = Tr(p(t)X n ), calculated with the true statistical operator, p(t). This is 
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achieved using Lagrange multipliers, which makes the relevant distribution to depend explicitly on the true statistical 
averages of the relevant variables, and thus on the true statistical operator pit) which is obtained by solving the 
modified Liouville equation |lj . This makes the whole scheme highly non- linear (and thus the term "self-consistency 
conditions"). 

When considering, within the MQCE, a mixed quantum-classical system consisting of quantum and classical 
degrees of freedom, one has to use a generalised expression for the Liouville operator which acts on classical variable 
A as the classical Poisson bracket 

where Q — {Qj} and P — {Pj} are the coordinates and momenta corresponding to the classical degrees of freedom, 
while it is the quantum commutator, ■h\) r , H], when it acts on a quantum operator Y associated with the quantum 
degrees of freedom. Here H is the system total Hamiltonian, depending on both types of variables at the same time. 
Consequently, the true and relevant statistical operators which may also depend on both classical and quantum 
variables, act as operators on quantum states and, at the same time, are functions of coordinate and momenta of the 
classical degrees of freedom, as in classical statistical mechanics. Since the quantum operators may not necessarily 
commute with the Hamiltonian, one has to use the symmetrised Poisson bracket when constructing the appropriate 
generalised Liouville operator |15 [ \W \ IT7 |, \E\ [24] . i.e. the Liouville operator in MQCE is the sum of the quantum 
and classical counterparts: 

it . . . = it, . . . + iL c . . . = 1 [. . . , H] + 1 ({. . . , H} - {H, . . .}) (3) 

i.e. it is formally constructed as a sum of the quantum and symmetrised classical Poisson brackets. It is readily seen 
that the Liouville operator defined in this way serves as a usual commutator when acting on quantum operators 
and is the classical Poisson bracket ([2]) when acting on classical variables. If a variable contains both classical and 
quantum components, the generalised operator ((3| is to be used. 

The statistical operator in MQCE is normalised to unity in the generalised sense, Tr (p) = 1, via the "total trace" 
defined as: 



Tr (...) = J tr(...)rfT (4) 

where the trace written with small letters corresponds to the usual quantum trace taken with respect to the quantum 
states associated with the quantum degrees of freedom, while integration corresponds to all coordinates and momenta 
of the classical phase space T as in ordinary classical statistical mechanics. Correspondingly, a statistical average of 
an arbitrary observable A, which may depend on classical degrees of freedom and, at the same time, is an operator 
in the quantum subspace, is defined in the generalised sense as 

(A)* — Tr (p(t)A) (5) 

One can also define the "flux" A operator (time derivative of A) associated with the variable A in the usual way as 
A = iLA. It is seen that all equations look identical to either pure classical and quantum cases, only the actual 
meaning of the Liouville operator is different. 

In the following, we shall limit ourselves with the Hamiltonian of the form 

H = T.^if + u ^) + H ^ ( i^) ( 6 ) 

3 3 

which corresponds to the electron-ion system in which ions of masses Mj are considered classically, while electrons 
quantum mechanically. The index j corresponds to a classical degree of freedom (i.e. each atom contributes three 
such degrees of freedom). Above, the first term gives the kinetic energy of the classical ions, their potential energy 
in an external field as well as the ion-ion interaction is provided by the second term. The last term forms the 
quantum Hamiltonian for the electrons (with coordinates q = {r^} and momenta p), 

H q (p,q;Q)=H e (p,q) + $(Q,q) (7) 

This Hamiltonian describes kinetic and interaction energies of the electrons (the first term) , while their interaction 
with the classical coordinates, Q, is described by the second term. Note that interaction between the two subsystems 
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in H depends only on their coordinates. Then the following expression for the classical part of the Liouville operator 
acting on the operator A is obtained: 



iL r A 



E 



Pj dA 
Mi <)(.), 



dA 
dP 3 



P + Pi 



dA 



(8) 



where the variable Pa — iLP^ 



iL c Pj 



OH 
dQj 



corresponds to the force acting on the classical coordinate Qj 



(degree of freedom j). 

When deriving kinetic equations (equations of motion for the statistical averages (A) 1 
proves to be indispensable: 



Tr 



A) = -Tr 



((ilA) P 



the following identity 



(9) 



where the generalised Liouville operator of Eq. ([3]) is used. Since this identity is linear with respect to the Liouville 
operator, it can be proven separately for each of the Poisson brackets. For the quantum bracket, iL q , it follows 
trivially from the cyclic invariance of the quantum trace [12J. To prove it for the classical Poisson brackets, one 
writes: 



Tv((iL c p) a) =^2 far 



El 

M, 



-tr 



dp 



A 



1 ( dp ■ , 



i 



:tr P 



8p_ 

J dP, 



A 



Using integration by parts with respect to Qj for the first term in the square brackets and the fact that the density 
operator should vanish at the boundaries of the phase space, we find that we can replace the trace tr (^-J^j-Aj with 

— tr^J^j-p). Similar method is applied to the other two terms: using integration by parts with respect to Pj, 
cyclic invariance of the quantum trace and the fact that ionic forces Pj do not depend on the momentum Pj , we 
obtain the following substitutions for the second and the third traces in the square brackets above: tr (^J^PjA 

— tr (pj§^pj and tr (^Pj-^-Aj — * — tr (j^prPjpj • This proves Eq. ([9j) for any quantum-classical operator A. 

Finally, we prove that for any general operator B(P,Q), acting on quantum states and depending on classical 
variables as well, the following identity is satisfied: 



Tr 



Ul c b) 







(10) 



provided that the operator B vanishes at the boundaries of the classical phase space. This is proven by using an 
explicit expression for the classical Liouville operator, Eq. |(8]). Indeed, consider the first part of it, containing 
the product of the classical momentum Pj and the derivative MM-. When taking the trace, the integration over 



Qj is performed immediately resulting in the difference B(P, Q) 



Qi — oo 



B(P,Q) 



which is zero due to our 



assumption concerning the operator B. Similarly, the other two terms of the classical Liouville operator, Eq. JH]), 
also result in the zero contribution due to integration over Pj and the fact that the force, Pj, does not depend on 
the momenta. 



3 Theory 



3.1 Relevant variables and distribution 

As we are primarily interested in this work with the equation of motion (EoM) for the classical variables which 
are much slower than their quantum counterparts, it is reasonable to sample over the fast degrees of freedom. 
In practice, this is achieved by choosing classical coordinates and momenta Q and P as the appropriate relevant 
variables. Correspondingly, the relevant distribution maximising the information entropy at the given temperature 
T and number of electrons N and subject to the self-consistency conditions, 



(Q J ->t eI = <Qi>*. ( P i)iel = ( P i) 



is (cf. [12] 



Prel (t) = — CXp { -f] 



H-^iVjPj + FjQj 



(11) 



(12) 
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where (3 — 1/kgT is the inverse temperature, H = H — p,N is the system Hamiltonian containing explicitly the 
chemical potential p, of electrons, Vj and Tj are the corresponding Lagrange multipliers and Z is the normalisation 
factor ensuring that Tr {p re i ) = 1 • Note that the sum over TV is incorporated into the definition of the "small trace". 
The relevant statistical operator depends on time only via the Lagrange multipliers (see below). 
At this point it is also convenient to introduce the statistical operator, 



Peq 



= Z- 1 exp(-pH g ) 



(13) 



for the quantum subsystem, Z eq — tr (exp (—(3H q )), where H q = H q — pN is the total electronic Hamiltonian ((7|)- It 
corresponds to the quantum equilibrium canonical statistical operator for the electrons when all classical variables 
are fixed (i.e. the classical subsystem is frozen). Then, the relevant distribution can also be written as a product 



of a "classical distribution function" 



f(P,Q,t) = ^exp{-0 



eq 
Z 



E 



Prei = Peqf(P,Q,t) 



\ - VjPj TjQj 



2M 



U{Q) 



J f(P,Q,t)dT = l 



(14) 



(15) 



and the canonical quantum equilibrium statistical operator, p eq . We shall also need the reduced distribution 
function, f(Q,t), which is obtained from the distribution function above after integrating over the momenta: 

f(Q,t) = J f(PQ,t)dP = ^3-e-P( u -^r j Q i ) ) J f(Q )t ) dQ = l (16) 

where Zq is the corresponding normalisation factor. Note that Z eq = Z eq (Q). The average with respect to the 
relevant distribution of any classical variable (depending only on classical degrees of freedom) is obtained as the 
average with respect to the distribution function /(P, Q,i); if the variable depends only on the classical coordinates, 
then the relevant distribution average is expressed as the average with respect to the Q-only distribution f(Q,t). 

The Lagrange multipliers are obtained from the self-consistency conditions (fTTIl . Due to explicit dependence of 
p re i on the ions momenta via f(P,Q,t), calculation of {Pj) re i is straightforward: it gives simply (Pj) rel = MjVj = 
{Pjf , and we obtain: 



V, 



(31 

Mi 



(17) 



Thus, Vj has the meaning of the average velocity of the degree of freedom j. 

Integrating with respect to all classical momenta in (Qj) rel , one finds that (Qj)* ei = (Qj) does only depend on 
the Lagrange multipliers {fj}. Inversely, this means that the Lagrange multipliers Tj only depend on the average 



coordinates |(Qj)*|. To obtain an explicit expression for Tj, we calculate (^Pj 
To this end, we first take the quantum trace of the identity (see, e.g. [12] ) 



rel 



§§r ) , using Eq. (L2 



dQj 



rel 



d 
dQj 



-fJH 



-a I , 
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dQj 



(18) 



giving 



\ tT ( JL e -PH ) = tr 



dH 



-0H 



\dQj~ J "\8Qj 
Therefore, using this and the integration by parts with respect to Qj in the 



(19) 



dH 
dQ 3 



one obtains: 



P, 



rel 



(20) 



Thus, the second Lagrange multiplier, Tj, has the meaning of the minus average force acting on the classical degree 
of freedom j. Using the explicit expression for the Hamiltonian, Eq. ([6]), the ionic force Pj can be broken down 
into conservative, 

,H (Xj) eq (21) 



dQj 
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and the stochastic, AXj — Xj — (Xj) eq , parts, Pj = F? + AXj, where Xj = — -§gp is the instantaneous force acting 
on the degree of freedom j due to the electronic (i.e. quantum) subsystem. Therefore, using Eq. I|14p . one obtains 
that 

Ti = ~ J f(P, Q, t) (P 3 ) eq dT = -J f( p , Q, t)F?dT = - (F?)l el (22) 

which demonstrates that the Lagrange multiplier Tj corresponds to the average of the conservative force. 

By definition, (AXj) eq = 0. However, using Eq. (HI), one can also check that the average of the stochastic part 
of the force with respect to the relevant distribution is also zero: 

(A^-)* ei = (23) 

The relevant distribution has a number of properties which are proven to be useful in our forthcoming analysis. 
Firstly, since p re i is equal to a product of a part, depending only on classical variables, and the quantum operator, 
e -l3H q ^ commu t es with the Hamiltonian, PL, Eq. J6]). Therefore, 

iL qPre i = and e lL " 1 p re i(t') = p re l{t') (24) 

Next, consider the relevant distribution average (Pjip(Q)) , where ip(Q) is some function of the classical 

v rel 



coordinates. Using the explicit expression for p re i, Eq. Ifl2j) . trace identity (fT9|) and integration by parts, one 
obtains: 

Note that Pj inside the angle brackets in the left hand side of this formula can be replaced with Pj since Pj = 
Fj + AXj, and the quantum equilibrium average of the stochastic force is equal to zero: 



{AX i>Y rel = J iP(Q)f(P,Q,t)tr(p eg AXj)dr ~ {AX 3 ) eq = 
due to Eq. lfl4|) . Therefore, by taking ip = Qi and ip = F[, we obtain the following two useful relationships: 

1 / B P c \ * 

(F? F ?Li ~ ( F ?Li <*?>'«! = — p {^-) rei (26) 

Since the left hand side in the second identity is symmetrical with respect to indices i and j, we also have the 
symmetry relation for the relevant average of the derivative of the conservative force: ^§^f~^) = 
Another useful expression is obtained by differentiating both sides of 

{Qi)* = {Q l )l el = J tr (QiPrei) dT = J Qif(P,Q,t)dT 

with respect to (Qj) ■ Recalling that only the Lagrange multipliers {Tj} in f(P,Q,t) (both in the exponential and 
in the Z) depend explicitly on (Qi) , we obtain: 

hij = y -^^dj (27) 

P jd{Qi) 3 V ' 

where C = (Q*Qj)* e j — (Qi)' {Qif ls a symmetric matrix. It follows from this that the derivative = q^q 3 )* 

is also symmetric. Because of Eq. l(22|) . the derivative of {Tj) , with respect to (Qi)* is also symmetric: 

ami, d( F ^ rel 

d(Qjf d(Q i ) t 
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To obtain a perturbative solution of the Liouville equation in Section I3.3[ some other expressions involving the 
relevant distribution are needed. By differentiating p re i of Eq. lfT2|) with respect to the classical coordinates and 
momenta, one gets: 



dP„ 



Mi 



dprel(t) 

dQ 3 







J J "3 

r,) +f: 

rel 



AXJi\h)d\ 



Prei 



(28) 
(29) 



where we have used Eqs. (fT?) . lfT8|) and lj22|) . and AXj(x) = e txn <i/ h /s.Xje~ lxHq / h is the stochastic force in the 
Heisenberg representation. 

Finally, one can also calculate derivatives of p re i with respect to the Lagrange multipliers, remembering that 
the "partition function" Z also depends on them: 



dVj 
dT, 



-■ (Pj - (Ptf) Prel(t) 
P (Qj ~ (Qj}*) Prel(t) 



(30) 
(31) 



3.2 Equations of motion for ions 

In order to derive EoM for ions, we should calculate the time derivatives of the exact averages (Pj) and (Qj) . 
Using the Liouville equation with broken time-reversal symmetry, Eq. ([1]), and identity |(9]), one obtains: 

!(<W'-*(^)-*(,(rt«.))-££ (32) 

where we have also used the fact that iLQ n = iL c Q n — P n /M n . Note that the term in the right hand side of the 
Liouville equation ([I]) does not contribute due to the self-consistency conditions (fTTjl . Similarly, using Eq. (J9j) , we 
get: 

j t (P n )* = Tr (p (iLP n )) = (Pn)* = (Pn)* + Tr (P n Ap) (33) 

where Ap = p — p re i is the difference between the exact and the relevant statistical operators. The obtained EoMs 
are similar to Newtonian ones of the ordinary Molecular Dynamics [H [2] since the right hand side of Eq. ([33)1 

corresponds to the actual force (Pnj acting on degree of freedom j. However, this force depends, in a rather 

non-trivial way, on the time evolution of the exact statistical operator p(t) which satisfies the Liouville equation 
([J). Note that p(t) is the statistical operator for the whole system, comprising both nuclei and electrons, and no 
attempt has been made to factorise p in any way here. 

The ion force (P n ) calculated with respect to the relevant distribution depends explicitly on time via the 

\ / rel 

Lagrange multipliers; the latter are some functions of the exact expectation values (Pj) and (Qj) as it has been 
discussed above. Similarly, the trace in the right hand side of Eq. ((33j) would depend on the Lagrange multipliers 
and thus on the same expectation values. One can see that the kinetic equations written above correspond to some 
non-linear differential equations for the observables (Pj) and (Qj) . To obtain these equations in the explicit form, 
we have to obtain an explicit expression for Ap(t) by solving the Liouville equation JT]). This will be done in the 
next subsection using a kind of a perturbation theory (cf. Refs. [14} [13] ) in which the square root of the relative 
mass of the electron and nucleus, y/m/M, is used as a small parameter. 

3.3 Perturbative solution of the Liouville equation: the first order 

Formally, the exact solution of the Liouville equation JT]) with respect to the Ap(t) can be written as [12] 

Ap(t) = - f dse" s e lsl (£. + if) p rel (r) (34) 
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where r = t + s. Here L is the combined Liouville operator, Eq. (J3j) , containing both classical and quantum parts. 
In order to apply the perturbation theory, we have to calculate the quantity + iLj p r ei(r). This calculation 
consists of several steps which will be outlined below. 

Firstly, iLp re i — iL c p re i due to Eq. (|24|) . The action of the classical Liouville operator J8]) on the relevant 
distribution is obtained using Eqs. (|28| and f29|) and is as follows: 



iL c p rel (r) 



^ Mi 



Pj / AXj(i\h)d\ 



Prel{r) 



(35) 



To calculate the time derivative of the relevant distribution, we note that it comes entirely from the Lagrange 
multipliers. The latter depend on time through the observables (Pj) and (Qj) which satisfy Eqs. f33|) and (|32j) . 
respectively. Also, we note that if Vj depends directly only on {Pj}*, the other Lagrange multiplier, Tj, depends on 



all average coordinates |(<3j)*|- Therefore, 



dpreljr) 

dr 



E 



/ dpreijr) 8Vj dpreiir) dT 
\ dVj dr dTj dr 



E M . ^ 



(Ps) r ) 



Pi 



Tr 



(PjAp(r))] P rel(r) 



{QjY){PiY Prel{r) 



(36) 



where the use have been made of Eqs. f30|) - ([33j) as well. 

Thus, the required quantity + iL \ p re i{r) is available now as a sum of two expressions, l|35p and (|36| . given 

above. These are to be acted with the exponential Liouville operator, e lsL (see Eq. lj34j) ). In turn, the Liouville 
operator consists of the quantum, iL q , and classical, iL c , parts. Similarly to the argument of Refs. [TSl [T4], we 
argue that the classical part of iL can be considered as being much smaller than its quantum counterpart. Indeed, 
on average, one can assume that the classical momentum Pj is of the order of Af 1 / 2 , where M is a characteristic 
mass of the ions. Then, according to Eq. ([8]), iL c ~ M~ x / 2 since the forces Pj depend only on ionic positions, 
not on their masses. Hence, the exponential operator, e zs ^ — e s [ lL i+ lL <=] ^ can De expanded in a power series with 
respect to the "small" Liouville operator iL c . This can be done systematically (and expressed via the time-ordered 
exponential operator, see, e.g., Chapter 6.1.1 in [23J); here we shall only need the first two terms: 



—ixsL c 



dx 



(37) 



In this Section we shall limit ourself with the very first term in this expansion, i.e. we replace the combined Liouville 
operator in Eq. (f34|) with its quantum part. Then, since the action of the quantum exponential Liouville operator 
on any quantum operator A is simply equal to its Hermitian conjugate, e lsLq A — A(s), we obtain in this order of 
the perturbation theory, using Eqs. i|35p and ((36j) : 

M*) = - J\ ds e" |E ^ [ (Fj ' r ^ + {P] {PjY) Tr ( Pj Ap{r) 



+Pj / AXj(i\h + s)dX 



P^ r ) ~ E ( p > - ^) r ) (^Xj(s)p rel (r) + Prel (r)AXj(s) 



y — (Pi) 



d(QiY 



(Qj - (Qj) r ) Prel(r) 



(38) 
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It is seen that this is actually an integral equation for Ap since it is present inside the integral (and the trace) as 
well. Therefore, it can be solved by iterations. To do this, we must analyse every term in the above expression with 
respect to the small parameter of our perturbation theory. It follows then that each term, apart from the term with 

the trace, is of the order of P/M ~ M" 1 / 2 (because of Eq. ([27|) . the derivative d g^q.y l should be considered as 
of the zero order with respect to M). Thus, in order to obtain Ap(t) in the first order, one can simply drop the 
trace term in the right hand side. We then immediately see that this particular approximation for Ap, which will 
be referred to in the following as Api, does not break down the normalisation of the density operator since, as it 
can easily be checked by direct calculation, Tr(Api) = 0. 

Hence, dropping the trace term above, substituting the resulting expression for Api into the kinetic equation 
([33]) . and recalling that the force P n = F° + AX n , we get two terms: Tr (i^Api(i)) = / F£tr (Api) dT and 
/ tr (Ap\{t)AX n ) dT. The quantum trace of Api, needed for the first term, is obtained directly from Eq. lj38j) by 
using the following identities: 

tr (/Vei(r)) = f(P,Q,r) 
tr(p rei (r)AX fe (s)) = f(RQ,r)tr(p eq AX k ) = 
where the cyclic invariance of the trace was used to obtain the second identity. Thus, we obtain: 



tr (Apr (t)) = - £ ds e" £ A J {Pj y _ <i^) - £ {Pj )> 



d{Q l ) r 



(Qi-{QiY)}f(P,Q,r) 



and, therefore, 



Tr (^Api(i)) = - ds e» £ JL (Ptf { ({F^)^ - (Ff)^ 
d(F c ) r \ 



This latter expression can be greatly simplified by virtue of Eqs. (|25| and (|26| : 



Tr(F„ c Api(i)) 



ds e £ 



dFJ 

9Qn 



d(F9) 



J I rel 



d(Q n ) r 



(PiY 



(39) 



Note that both derivatives inside the square brackets are symmetric with respect to the permutation of their indices. 

We shall now turn to the second term, Tr (ApiAX„), arising in the right hand side of the kinetic equation, and 
substitute Api there. Noting that the trace Tr \FjAX n p re i) is equal to zero due to tr (AX n p eq ) — 0, and the fact 
that 

(40) 



(41) 



PjPreli^dPj = {PjY j Prel{r)dPj 

which follows from the explicit dependence of p re i on the classical momenta and Eq. ifTTj) . we obtain: 





Tv(A Pl AX„ 



[ dse™J2 

J — OO 



(P 3 ) r (X n ,X 3 (s)) 



rel 



where we have introduced the non-equilibrium correlation function of the fluctuation of the ionic force (cf. [T4l[T3] ): 



(X ni X k (s)) 



rel 



(AX n AX k (ih\ + s)) d\= dQf(Q,r) 



rel 



(AX n AX k (ihX + s] 



dX 



eq 



(42) 



The last passage in the above formula is due to the fact that the relevant distribution average under the A— integral 
depends only on the Q variables since the P integration can be performed directly. Combining Eqs. (|39|) and f4Tj) . 
we finally obtain the kinetic equation for the ionic momenta in the following form: 



dt 



nl rel 



E 



ds e e 



dFl 
dQ r , 



d ( F iY 

V J I rt 

d(QnY 



P{X n ,Xi{a)l 



Ah 



(43) 
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This is the desired equation of motion for ions. In the right hand side it contains the total force acting on ion n 
due to other ions. Their interaction and energy exchange with all the electrons are also completely accounted for. 
We also note that we have not made any assumptions as to whether electronic subsystem is in its ground state, 
it is in general a weighted sum of the ground and excited electronic states (see also below). In other words, this 
description goes beyond the adiabatic approximation. 

The first term in the right hand side of Eq. |43|) gives the conservative force corresponding to ionic positions at 
time t: 

(K)ll = I F c n {Q)f{Q,t)dQ (44) 



The second term in the right hand side of Eq. i[43|) gives a correction arising due to fluctuation of the ionic force. 
Similarly to the friction force acting on a Brownian particle immersed in a liquid, this force appears to be linear with 
the ions momenta. Indeed, the correlation function and derivatives inside the square brackets do not depend on the 
expectation values (Pj) since the P integration in those terms can be performed explicitly and all the Lagrange 
multipliers {Vj} disappear exactly. Moreover, the derived friction is non-Markovian, i.e. includes memory effects. 

Thus, we conclude, the rigorous non-equilibrium statistical mechanical treatment of a system composed of ions 
and electrons results in Newton-like equations of motion for average ionic momenta that additionally contain friction 
forces due to energy exchange with the electronic subsystem maintained at the given temperature T. 

The obtained equations are very complicated because of the relevant distribution used in the right hand side 
which depends on the observables (Pj)* and (Qj) in a rather complicated way. In the next subsection a reasonable 
approximation will be offered which results in a significant simplification of these equations. 

3.4 Saddle-point approximation 

Let us consider a relevant distribution average of some function, C(Q), depending only on classical coordinates: 

(C(Q)) t rd = J C(Q)f(Q,t)dQ 

where the Q-distribution is given explicitly by Eq. lfT6|) . Let {ip m (q,Q)', m = 0,1,2,...} is the complete set of 
electronic wavefunctions, depending parametrically on the positions of ions Q. The wavefunctions ipm are the 
eigenvectors of the electronic Hamiltonian, i.e. H q il> m = (^m — A*) V'm- Then, it is easily seen that the "partition 
function" of the Q— distribution, f{Q,t), can be written as a sum: 



Z Q = J e-W» I 1 + J2 g— /3Ae m (Q) I dQ (45) 

where 

R(Q) = e (Q) — fJ> + U(Q) + £ (ff)*^ Qj (46) 

j 

and Ae m (Q) = e m (Q) — £q(Q) are exact electronic excitation energies for the given geometry of the nuclei, Q. We 
assume hereafter that the ground state is non-degenerate for any geometry, and thus all the excitation energies are 
strictly positive. Moreover, we assume that for any geometry Q there is a gap between the ground and the first 
excited states, and the ground state energy in the external field, £o(Q) + U(Q), has a minimum at some geometry 
Qo- Of course, the minimum will be affected by the last term in Eq. {46]), however, we assume that this term does 
not change significantly the potential energy surface of the ground state. Therefore, the function 3?(<5) will still 
have a minimum at some geometry Q* (the subscript reflects the fact that, because of the last term in Eq. lj46j) . the 
minimum geometry Q l will depend on time), and thus can be expanded in a series with respect to the difference 
Q-Q': 

W) - + 2 E {9Q-9Q-) Qt («« - ) - «5) + • • • 



where the matrix 



( d 2 n \ 

\dQ i dQ j ) ( 



of second derivatives is positively defined (since Q l is the minimum). Hence, the 



function e~^^) will be highly peaked around Q* , whereas the function in the round brackets in Eq. f45|) can 
be assumed to be a rather slowly changing with Q and can thus be taken away from the Q-integral with all ionic 
positions calculated at Q = Q*. A simple calculation in the spirit of the well-known saddle-point approximation 
will then show that the distribution function f(Q,t) effectively serves as a Delta function S (Q — Q*) giving for the 
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average of any slowly changing function of ionic positions, C(<9); the following simple result: {C(Q))rei — C(Q*)- 
Further, if we consider specifically C(Q) = Qj, then we obtain that (Qj) 1 — (Qj)* e ; should be replaced with Q*. 
Thus, we conclude that in a consistent application of the saddle point approximation, one replaces the averages 
(C(Q))rei) calculated with respect to the Q— distribution, with the corresponding functions £ calculated at exact 
average ionic positions at time t, namely: 

(47) 



In particular, this result can be applied to the non-equilibrium correlation function of Eq. (|42|) resulting in a 
much simpler expression for it: 



(X n ,X k (s))l el 



AX n AX k (ih\ + s] 



dX 



= (X n ,X k (s)) { e f 



(48) 



-I Q={Q) t 



Thus, the correlation function depends directly on a single time s; however, the statistical average is to be calculated 
over the equilibrium distribution p eq corresponding to ions fixed in their exact positions (Q)* at another time t. 

Using the same approach, one can also verify that the term in the square brackets in the right hand side of the 
equation of motion (|43| can be dropped. Indeed, since the conservative forces F? depend entirely on ionic positions, 
we can write: 



dF[ 



dF[ 
dQ n 



Q=(QY 



d(Q n ) r 



and, at the same time, 



d ( F !Y rel 9F9({QY) 



d(Q n y 



d(QnY 



Therefore, within the saddle-point approximation, the difference of derivatives in the square brackets in the equation 
of motion f43|) is equal to zero. Following the same arguments and replacing the conservative force {F^)^ in the 

equation of motion with F£ ^(Q)*^, we obtain: 



dt 



F. 



n (<<?>') 



^ Mi 



dse^ (X„,^-( S ))<f (P 3 ) r 



(49) 



J J-oo 



which is the final result. 

We see that, if not for the friction term, the equations of motion would have corresponded exactly to the 
Newton's equations of motion for ions: in the left hand side we have the time derivative of ion n momentum, while 
in the right hand side - the total statistically averaged force acting on this ion at the given temperature: 



F, 



n ((Q)') 



dH 

dQ n 



Q=(QV 



dU 

9Qn 



(ipr, 



dH q 

dQn 



\lpr 



q=(qY 



dU 
dQn 



-/3(e m -M) 



de r , 



dQ r 



(50) 



Q=(Qy 



fi)). Here, — fn 2 - is the force acting on ion n (due to electrons and nuclei) when the 



i{Q), corresponding 



where Z eq = J2 m ex P i~P ( £ m rjj- -g^ 

electronic subsystem is in electronic state m (i.e. on the adiabatic potential energy surface, e r 
to electronic state m). 

Note that atomic positions, (Q)*, correspond exactly to the averaged ionic momenta, (P) , see Eq. (|32j) . The 
obtained equations of motion are more general than those of ordinary Molecular Dynamics. Indeed, in standard 
MD the forces do not depend on temperature and are calculated as in Eq. f50|) from the conservative part taking 
only the electronic ground state into account, 



F c ~ - 



dU 
dQ n 



<V>o| 



dH q 

dQn 



d(U + e ) 

dQn 



(the Car-Parinello ground state ab initio MD simulations [4] being an obvious example). Therefore our equations 
may serve as a justification of MD simulations which go beyond the Born-Oppenheimer approximation (see, e.g. 
[TT]). where, when calculating the force acting on an ion, it is assumed explicitly that the electronic subsystem may 



11 



occupy both ground and excited electronic states. In fact, in complete agreement with the principles of quantum 
statistical mechanics, we show that there is a certain probability for the electronic subsystem to occupy every state 
at the same time, which, one must admit, is somewhat different from most of the non-adiabatic computational 
techniques (6j in which it is usually assumed that only one state can be occupied at every single time step. 

We also observe that a consistent non-equilibrium treatment results in additional terms in the equations of 
motion which are proportional to the ions momenta and thus have the meaning of friction, related to the energy 
exchange between the ions and the electrons; the latter serving as a "thermostat" held at a given temperature. 
Thus, our rigorous treatment justifies the usage of "electronic friction" terms in MD simulations [IT] and explains 
their physical origin. 

3.5 Non-equilibrium correlation functions 

Since the operator of the atomic force, Xi — —-§§- — J2k=i x i( r k), is a derivative of the electron-phonon interaction 



energy, Q) = X^fcLi <K r fc! Q)> it is a one-particle operator, 

ab 

(where c\ and Cf, are creation and annihilation operators in some basis set of spin-orbitals) and X % ah — (a \ Xi(r) \b). 
Here Xi(v) — — is the force on atom i due to a single electron at r, 4>(r,Q) being the interaction energy of 

this electron with all nuclei. Hence in general, the correlation function l|48p is a two-particle equilibrium statistical 
average, containing four c-operators, and thus cannot be calculated exactly in the general case. 

The calculation is straightforward if the Hamiltonian H q — J2ab habctcb is a one-particle operator (e.g. in the 
Hartree-Fock approximation). Indeed, in this case one can diagonalise the Hamiltonian, 

a 

where 

da = e laCa, 4 = E e o a c\ (51) 

a a 

and £o- and e CT = \\e aa \\ are the eigenvalues and the eigenvectors of the matrix h = \\h a b\\, and therefore express the 
operators Xi and Xj via the operators d\ and d a , e.g. 

X i = ]T a&dUr', 41 = E Kb^a'b (52) 

era' ab 

Since the calculation of the operators d} a {t) and d a (t) in the Heisenberg representation is simple, we obtain: 

(X,, X^t))^ = ]T x%x%n a (1 - tv) X ~ k) e^-'-W*/* (53) 

aa' 

where n a = {e^a-v) + ly 1 and x ( E j = E -i (i _ e -/3B) if E ^ and x(0) = (3, Note that the superscript (Q)* 
to the correlation function has been omitted to simplify the notations. 

When the electron-electron interaction is accounted for explicitly, one has to develop more powerful methods. 
First we note that the correlation function, Eq. (|48j). may be considered as a particular case of a more general 
correlation function defined for any two quantum operators A and B as follows: 



{A(t 1 ),B(t 2 )) e = [ dAtr[!(ti)B(t2+iA%J = /' dx/A^B^+iXh)) (54) 
Jo 1 J Jo x 1 e i 

This correlation function obeys some simple symmetry properties: 

(A(ti), B(t 2 )) eq = (A, B(t 2 - h)) eq = (A(t 2 ), Bitt))^ (55) 

which follow from the cyclic invariance of the trace (change of variables A — > Ai = /3 — A is also necessary to obtain 
the last equality). 
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An explicit expression for the correlation function can be obtained using the complete set, {ipn}, of the eigen- 
vectors of the electronic Hamiltonian (i.e. H q ip n — E n tp n ): 

(A, B(t)) eq = J2 PnX(E m - E n ) (ip n \A \if> m ) B \$ n ) e'^-emJt/R (56) 

nm 

where p n = Z~ q 1 e^ l3 ^ En ~^ N \ It is also convenient to introduce the spectral function (matrix) as the Fourier 
transform of the correlation function: 

(A,B(t)) eq = -J^ t J AB (u;)doj (57) 

where 

Jab{oj) = 2it%y^ j p n x{E m - E n ) {ip n \ A \ip m ) (ip m \ B \tp n ) 6 (e m - £„ - fuS) (58) 

nm 

which satisfy the following symmetry properties: 

Jab{u) = Jba(-u) = Jba{u)* (59) 

Unfortunately, it appears impossible to develop systematically a perturbation theory for the direct calculation of 
the correlation function (|54| . It has been found, however, that such a perturbation theory exists for the equilibrium 

statistical average inside the A— integral, i.e for the (A(ti)B(t 2 + i\h)\ . Once the statistical average is calculated, 

\ / eq 

the correlation function follows immediately. The method is based on a relationship between the correlation function 
above and the Matsubara Green's functions (e.g. |22l I23|); this is revealed in Appendix. Therefore, one can use 
the perturbation expansion (and thus the powerful diagrammatic techniques) to account for the electron-electron 
interactions in calculating the correlation function appearing in Eq. lj48j) . 



3.6 Perturbative solution of the Liouville equation: second order 

The formulae developed in the previous Sections correspond to the first order approximation as we kept only terms 
of the order of M -1 / 2 . This treatment can be systematically extended to higher orders. Unfortunately, this results 
in very cumbersome expressions even in the second order. Therefore, in this Section we shall simply outline the 
main idea of how the extension to higher orders can be done. 

We start from Eq. 11341) and replace the exponential operator there with its expansion given by Eq. lj37j) . After 

that, we recall that the operator + iLj p re z(r) contains the trace of Ap as well (see Eq. (J36J) ) and thus this 

dependence must also be included in developing the perturbative expansion. For instance, one gets for the next 
order term: 



Ap 2 (t) = - J ^ ds e" £ A ( P . _ {Pj y) prel {r ) Tr (P 3 A Pl (r)] 

dse £S / dxe lxsZ " (isL c ) e - ia!flZ «e _ " £ «A (60) 
> Jo * ' 



where A is obtained by removing the trace term in the operator + iLj p re i{r) (see Eqs. I|35p and f36|) ) and 

Ap\ is the first order term given by Eq. |38j) . 

First of all, it can easily be shown that Tr(Ap2) = 0, i.e. this correction is consistent as well with the correct 
normalisation of the statistical operator. Indeed, the trace of the first term in Eq. f60|) vanishes due to integration 
over Pj in the classical part of the trace and Eq. (J40j) . The second term also does not contribute to the trace of 
Ap2 due to the cyclic invariance of the trace and the operator identity (flC)|) (the operator A is proportional to the 
relevant distribution and thus vanishes at the boundaries of the classical phase space) . 

It can also be shown that any order correction to the statistical operator, Ap n , has zero trace and hence does 
not break down the normalisation of the statistical operator. 

Using the explicit expression for Ap 2 given above, one can calculate its contribution to the force, Tr (P n Ap2(t)j , 

in the right-hand side of the equations of motion ([33)1 . The contribution is very cumbersome (and is only due to the 
second term in Eq. l(60|) ) and won't be reproduced here. We only note that it is proportional to the square of the 
atomic momenta, i.e. it contains terms proportional to (Pj) r (PkY ■ Higher order terms contain more products of 
the atomic momenta. Also, much more complicated correlation functions appear as the kernels of the time integrals 
in the contribution to the force. 
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4 Conclusions 



In this paper we have considered, using intrinsically non-equilibrium statistical mechanical theory, a system of fast 
(electrons) and slow (atoms) particles which interact with each other and may interchange their energy. The system 
is enclosed in a thermal bath kept at a constant temperature. In our treatment, the combined quantum-classical 
consideration was used: the slow degrees of freedom (atomic coordinates) were treated classically, while the fast 
variables (electrons) quantum mechanically. No assumption was made concerning the structure of the statistical 
operator; in particular, it is not in any way factorised. In addition, the electronic subsystem was treated exactly 
with complete inclusion of electron-electron interaction. 

We show, by assuming that the classical degrees of freedom are much "heavier" than the quantum ones, that 
equations of motion for the former (i.e for atoms) contain the conservative and friction forces. The conservative 
forces are statistically averaged over the electronic states. The friction force, which is strictly proportional to the 
atoms momenta, is expressed via the correlation function of the fluctuating force with which electrons act on the 
atoms (i.e. due to the fluctuation of the electron-phonon interaction). The correlation function can be expressed via 
two-electron Matsubara Green's function and thus calculated using the well-developed perturabtion diagrammatic 
techniques. 

The theory presented here gives a solid foundation for a number of intuitive theories based on MD simulations 
which go beyond the Born-Oppenheimer approximation (see, e.g. [HI [6]). We also justify the usage of "electronic 
friction" terms in MD simulations [TT] and explain their physical origin. 

In our method it was assumed that on average the electronic subsystem is in thermodynamic equilibrium. For 
instance, the theory developed here is applicable to ordinary Molecular Dynamics when atoms move along classical 
trajectories. However, the theory can also be applied in other cases in which the electronic subsystem is in a 
steady state, i.e. which, on average (more precisely, over the timescale associated with atomic motion), is not time 
dependent (is stationary), for instance, stationary electronic or heat conductance in an atomic wire [HE]. 

A number of avenues exist in developing our theory further. Firstly, as has been mentioned above, we have 
assumed in our treatment that electrons quickly reach thermodynamic equilibrium during the motion of atoms. 
In some cases this assumption will not be valid, e.g. when considering non-elastic effects during a nonstationary 
conductance in a system. In these and similar cases one has to include additional relevant variables into the 
consideration, e.g. time dependent electronic density matrix as in the kinetic theories [12] . This approach would 
result in an additional kinetic equation for the density to be solved simultaneously with the atomic motion considered 
here. Secondly, another possible extention is concerned with considering nuclei quantum mechanically as well. This, 
however, is much more difficult and may be done by e.g. starting from the method developed in Ref. [9]. 
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Appendix 

In this Appendix we shall relate the correlation function lj48j) with the thermodynamic Matsubara Green's function 
(e.g. [221 [23]). We shall start by defining a "complex time" Green's function for two arbitrary operators A and B 
as follows: 

Gab{xi,x 2 ) — -tr p eq T x (A(xi)B{x 2 )) (61) 

where the quantum operators A and B appearing in the Green's function are written in the special representation 
defined with respect to a single real parameter x. This specific representation which will be designated in what 
follows with the bar over the operator symbol, is given by 

A(x) = e iT(x)H q /h Ae - t T( x )H q /h ( 62 ) 

where the complex "time" t(x) = (x is introduced which contains the complex prefactor ( = — (t + i\h) / A and 
changes linearly with x, so that r(0) = and r(— A) =t+ i\H. The operator T x performs chronological ordering of 
the operators A(x\) and B(x%), so that x increases from right to left: 

fx (A(x 1 )B(x 2 )) = 0(x 1 -x 2 )A(x 1 )B(x 2 ) + r ] 9(x 2 -x 1 )B(x 2 )A{x 1 ) (63) 
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where rj = ±1 corresponds to the sign acquired when changing the order of the operators A(x\) and B(x 2 ); 9(x) is 
the Heviside step function. 

Using the cyclic invariance of the trace, it is easy to see that the Green's function (|6lj) actualy depends only on 
the difference x — x\ — x 2 : 



G AB ( Xl ,x 2 ) = G AB (x) = ~9 (x) (AB{~x)) - r,6 (-x) (B(-x)A) 



<■<] 



The correlation function we would like to calculate, (AB(t + iXH)) with < A < (3, is equal directly to the 

\ / eq 

minus Green's function G AB (0,—\) = Gab (A). Moreover, using eigenvectors and eigenvalues of the Hamiltonian, 
H q 4>n = (e« — fiN) ip n , one obtains: 

(AB(-x)) eq = (i>n\Aty m ) {^ m \B\^ n ) e^™-^ xt / Xh e- x ^-^\ x > (64) 

nm 

(B(-x)A) eq =Y,Pn «>n| A \^ m ) <^ m | B \1> n ) e <(=»-««)«*/Aft e -(x+«(c„-«„) j x < ( 65 ) 

nm 

where p n = Z^e - ^^" - ^^. It is easily verified (by splitting each double sum into two contributions with positive 
energy differences each) that the first correlation function, Eq. ((64J) , converges for x < (3 (and x > 0), while the 
second one for x > —(3 (and x < 0). Therefore, x is limited to the interval —j3 < x < (3 . 
Note also that the above correlation functions can also be written as Fourier integrals: 



e iu>xt/\ e -x U h X>Q 



OO 

G AB {x) = -{AB{-x)) eq = J ^Iab(lo) 

— OC 
OC 

G AB {x) = - V (B{-x)A) eq = i 1 J ^-I AB (cj)e^ t / x e-^ + ^ n , x<0 

where the spectral function 

Iab{u)) = -27T^p„ (?p n \ A \ip m ) (ip m \B \ip n ) 6 . 



- LO 



So far, the Green's function introduced above in Eq. f61j) has been shown to possess properties very similar or 
even identical to those of the Matsubara Green's function [22, 23]. To strengthen this analogy, one can also expand 
the Green's function into a Fourier series in the interval —(3 < x < (3 or notice that the two Green's functions, for 
— j3 < x < and < x < (3, make a jump at x = 0: 



oo 

lim o [G AB (5)-G AB (-5)]= / — I AB {u){l-r)e 



However, this analogy is not complete; for instance, one can see from the integral representations of the Green's 
function given above that G AB (x) for x < is not related to the G AB (x + (3) although they both share the same 
spectral function I ab (lu). 

At this point we are quite prepared to find the relationship between the "complex time" and the Matsubara 
Green's functions. To this end, we shall first introduce the appropriate "interaction representation" for the operators: 

where Ho is the Hamiltonian of non-interacting electrons, i.e. H q = Ho + H' with H' being the electron-electron 
interaction. Then, the product of two operators in the Green's function f61j) can be written as 

A(x 1 )B{x 2 ) = U(0,x 1 )A I {x 1 )U(x ll x 2 )B I (x 2 )U(x 2 ,0) (66) 

where 

U(xi,X 2 ) — e iT ( x i) H o/tl e -i{T{x 1 )-T(x2))H q /h e -iT(x 2 )'Ho/h ^qy) 
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is the evolution operator in our interaction representation. The evolution operator defined above satisfies the 
differential equation (k = $ = 1 — 4v) 



dU(x 1 ,x 2 ) 
dxi 



-kH i (xi)U(x 1 ,x 2 ) 



which can be converted into an integral equation and then solved iteratively, giving as a solution: 



U(xx, x 2 ) = T x exp ( — k / Hj(x)dx 

J X2 



(68) 



Thus, the evolution operator can be expressed via the T^-exponent. Note that k actually depends on the particular 
values of t and A used to define the "complex times" in the Green 's function arguments. 

The fact that the evolution operator can be expressed via the chronologically ordered exponent allows for the 
considerable simplification of the chronological product of the two operators in the Green's function. Indeed, the 
product in Eq. f61j) can be rewritten as 

f x (A(xi)B(x 2 )) = f x (C/(0, x 1 )A I (x 1 )U(x 1 ,x 2 )B I (x 2 )U{x 2 , 0)) 

= f x (C7(0, x 1 )U{x ll x 2 )U(x2,0)A I (x 1 )B I (x 2 )) 

= f x (U(0, O^ix^Bjfa)) = f x (A I (x 1 )B I (x 2 )) (69) 

where use has been made of the obvious properties of the evolution operator: U(x,x) — 1 and U(xi, x 2 )U(x 2 , X3) = 
U(x\,X3). Thus, we see that the electron-electron interaction can be actually eliminated entirely from the chrono- 
logical product of the operators in the Green's function. 

Although the method developed below is valid for any operators A and B, we shall consider the particular 
case needed here when the operators are of the binary form given by Eq. (|52| (note also that r\ = 1 in this 
case). To proceed, we recognise that the creation and annihilation operators, S a and d a (in the representation that 
diagonalises the Hamiltonian Ho, see Section l33| . have a very simple form both in the standard thermodynamic 
(imaginary time) and our (complex time) representations. To simplify the notations, we shall use from now on in 
this Appendix a wavy line above operators for the usual thermodynamic interaction representation of the operators, 
i.e. d(x) = e xUa Ce- xUo . Then, one has: 

d aI {x) = d a e- lT ^°' r \ dl^x) = dU" (x)(Jh 

d aI (x) = e xHo d a e~ xn ° = d a e~ x ^, $ aI (x) = d) a e x ^ 

Therefore, any binary operator, A = J2aa' A aa >d'ld a > , when written in the "complex time" interaction representation, 
Ai(x), can easily be expressed as another operator A' written in the ordinary thermodynamic representation, i.e. 



At{x) =J2 A ™' d l( x ) d °'( x ) =y^K<7>{x)dt{x)d (T ,[x) =A'j{x) 



with the new matrix of the coefficients A'(x) = ||e^ K i"''>A tTtT i || which depends explicitly on x. This simple 

result allows us to rewrite the chronological product of the operators of Eq. I|69p as 



T, 



(A( Xl )B(x 2 )) = f x (A I (x 1 )B I (x 2 )) = f x (A'^xJB'^) 



(70) 



In the final expression above the product of two operators appears exactly as in the thermodynamic (Matsubara) 
Green's function. To finish the transformation, we should introduce the evolution operator in the usual thermody- 
namic representation: 



U(xi,x 2 ) 



e x 1 Ho e -(x 1 - a:2 )H ae -x 2 Ho =7; exp 



H' I (x)dx 



(71) 



that satisfies the properties U(xi, x 2 )U{x 2y , X3) = U(xi,x^) and U{x,x) = 1. The evolution operator enters the 
equilibrium statistical operator, p eq [221123] : 



Peq eq " 



-P{Ho+H>) = 



T x exp 



Pa 



tfH>(x)dx 



T x exp 



IoH'j(x)dx 



o=^PoU((3,0) 



(72) 
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where Z eq = Z§ (U ((3, 0)y and the brackets (...) = tr [po . . .) correspond to the statistical average with respect to 

pa = Z^e-W" with Z a = tr (e"^ ). 

The following steps depend on the particular values of the arguments x\ and x 2 . However, since the Green's 
function depends only on their difference, x = x\ — X2, which lies between -(3 and /?, it is sufficient to consider only 
negative values of x\ and x% . Thus, combining Eqs. (|70| and f72|) . we can write for the product of the operators in 
the Green's function i|6ip : 

p eq f x (A(x 1 )B(x 2 )) = ^ Po U(P,0)T x [^(zijff^afc) 

The evolution operator U(f3, 0) contains the sum of ordered products of operators Hj(x) whose arguments x lie 
between zero and (3, i.e. are all positive. Since, by our assumption, both X\ and x 2 are negative, the above formula 
can be transformed into: 

Z,: - \u(p 1 0)A' I (x 1 )B' I (x 2 ) 



-poT x 

^eq 

which results in the following final expression for the Green's function: 

o 



G A b{xi,x 2 ) = - 



T x U(p 1 0)A' I (x 1 )B' I (x 2 ) 



U(J3,Q) 



T x exp 



T^exp 



-f P H'(x)dx 



(73) 



(74) 



which is nothing but the Matsubara Green's function, Ga'B'(xi,x 2 ). The latter is defined with respect to the 
operators A' and B' which are obtained from the original operators A and B by using the primed matrices of 
coefficients as explained above. 

Thus, there is a direct connection between the Green's function lj6lj) and the appropriate Matsubara Green's 
function. Since there is a well-known diagrammatic perturbation technique developed for the latter with the 
denominator cancelling out exactly as a prefactor to connected diagrams in the nominator |22[ [23] , 



Gab(xi,x 2 ) 







/fa exp 


— \ H'j{x)dx 




Jo 



A' I (xi)B' J ( k X2) 



(75) 



where the subscript "c" indicates explicitly that only connected diagrams are to be retained, this method can be 
directly used to obtain corrections beyond the one-electron approximation. The latter was employed in Section l3~5l 
to derive the formula for the correlation function. In particular, in the zero order (when the exponential operator 
above is replaced by unity) , the same expression is obtained for the correlation function as in Section 13.51 Notice 
that the direct application of Eq. l(75|) results in an expression containing an additional term with the product of 
averages (A) {B} ; this term did not appear in Section [3751 since the correlation funciton considered there contained 
already the difference operators A A and AB. 
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